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I.  INTRODUCTION 


The  present  report  represents  the  latest  portion  of  an  ongoing  effort 
aimed  at  the  development  of  a  Navier-Stokes  analysis  for  the  general  steady  or 
unsteady  cascade  flow  field  problem.  Although  the  two-dimensional  cascade 
analysis  represents  a  simplified  version  of  the  actual  three-dimensional  flow 
field  which  includes  end  wall  effects,  the  two-dimensional  problem  gives 
significant  insight  into  the  cascade  flow  field  and  obviously  is  a  necessary 
first  step  in  developing  a  three-dimensional  analysis.  Hence,  cascade  analyses 
of  various  types  have  been  a  subject  of  high  interest  in  recent  years.  Among 
the  analyses  being  pursued  are  inviscid  analyses  (e.g.  Refs.  1  and  2),  inviscid 
analyses  with  boundary  layer  corrections  (e.g.  Ref.  3)  and  Navier-Stokes 
analyses  (e.g.  Ref.  4  and  5).  Each  of  these  approaches  are  viable  under 
certain  circumstances.  For  example,  inviscid  analyses  can  give  good 
predictions  of  the  blade  pressure  distribution  for  conditions  where  the  effect 
of  viscous  phenomena  upon  the  blade  pressure  distribution  remains  small. 
However,  inviscid  analyses  require  some  method  of  assuming  airfoil  circulation 
to  obtain  a  unique  flow  solution.  In  general  this  specification  is 
straight-forward  for  sharp  trailing  edged  blades  in  steady  flow  where  the 
Kutta-Joukowski  condition  serves  to  specify  circulation.  However,  in  cases 
where  the  trailing  edge  is  rounded  or  the  flow  is  unsteady  or  trailing  edge 
separation  occurs  specification  of  a  proper  condition  is  not  clear.  Finally, 


1.  Farrel ,  C.  and  Adamczyk,  J.:  Full  Potential  Solution  of  Transonic 
Quasi-3-D  Flow  Through  a  Cascade  Using  Artificial  Compressibility.  ASME 
Paper  81-GT-70,  1981. 

2.  Casper,  J.R.,  Hobbs,  D.E.  and  Davis,  R.L.:  The  Calculation  of  Two- 
Dimensional  Compressible  Potential  Flow  in  Cascades  Using  Finite  Area 
Techniques.  AIAA  Paper  79-0077,  1977. 

3.  Hansen,  E.C.,  Serovy,  G.K.  and  Sockol,  P.M.:  Axial  Flow  Compressor  Turning 
Angle  and  Loss  by  Inviscid-Viscous  Interaction  Blade-to-Blade  Computation, 
Journal  of  Engineering  for  Power,  Vol.  102,  1980. 

4.  Shamroth,  S.J.  and  McDonald,  H.:  An  Assessment  of  an  Ensemble-Averaged 
Navier-Stokes  Calculation  Procedure  for  Cascade  Flow  Fields.  Scientific 
Research  Associates  Report  R82-920011-F,  1982. 

5.  Shamroth,  S.J.,  McDonald,  H.  and  Briley,  W.R.:  Prediction  of  Cascade  Flow 
Fields  Using  the  Averaged  Navier-Stokes  Equations.  To  be  published  in  ASME 
Journal  of  Engineering  for  Power. 
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inviscid  methods  obviously  cannot  give  either  heat  transfer  or  viscous  loss 
estimates.  Nevertheless,  despite  these  inherent  limitations,  inviscid  analyses 
are  valuable  for  the  prediction  of  blade  pressure  distributions  for  flows 
having  little  viscous  displacement  effect  and  a  clearly  applicable  Kutta 
condition. 

Some  limitations  of  inviscid  analyses  can  be  relieved  by  combining  an 
inviscid  flow  calculation  procedure  for  the  pressure  distribution  with  a 
viscous  flow  boundary  layer  development  in  either  a  strong  interaction  or  a 
weak  interaction  mode.  A  recent  example  of  a  combined  viscous-inviscid 
procedure  is  the  work  of  Hansen,  Serovy  and  Sockol  (Ref.  3).  In  cases  where 
the  viscous  displacement  effect  has  an  insignificant  or  only  small  effect  on 
the  actual  blade  pressure  distribution,  an  inviscid  calculation  can  be  made  to 
obtain  the  pressure  distribution  and  a  boundary  layer  calculation  then  made  to 
obtain  heat  transfer  and  loss  effects.  However,  in  many  cases  the  viscous 
displacement  effect  may  significantly  alter  the  pressure  distribution.  Such 
cases  are  found  when  boundary  layers  become  thick  or  separate  or  in  transonic 
flow  where  the  local  pressure  distribution  and  shock  location  become  very 
sensitive  to  small  change  in  the  effective  passage  area.  In  these  cases  a 
strong  interaction  solution  is  required  to  account  for  the  mutual  effects  of 
the  viscous  boundary  layer  and  the  nominally  inviscid  core  flow. 

A  strong  interaction  analysis  may  take  the  form  of  either  a  forward 
marching  procedure  or  a  global  iteration.  For  regions  where  the  outer 
nominally  inviscid  flow  is  supersonic  (and  thus  described  by  hyperbolic 
equations)  a  solution  can  be  spatially  forward  marched  in  the  nominally 
streamwise  direction  with  the  inviscid  and  viscous  regions  coupled  on  a 
station-by-station  basis.  The  chief  difficulty  with  this  approach  is  the  stiff 
nature  of  the  coupled  sets  of  equations  which  is  manifested  in  the  appearance 
of  physically  unrealistic  branching  solutions.  In  regions  where  the  outer  flow 
is  subsonic,  the  outer  flow  equations  are  elliptic  in  nature  and  in  these 
regions  forward  marching  in  the  streamwise  spatial  direction  is  not  possible 
and  consequently  a  series  of  viscous  and  inviscid  calculations  must  be 
performed  in  which  each  corrects  the  other  in  a  global  manner.  Problems  with 
interaction  solutions  become  particularly  severe  in  transonic  flows  where  both 
subsonic  and  supersonic  nominally  inviscid  regions  are  present  and  where  small 
viscous  displacement  effects  may  have  a  major  effect  on  the  blade  pressure 
distribution  and  shock  location.  A  final  difficulty  with  the  interactive 
approach  occurs  when  boundary  layer  separation  appears.  Here  with  an  imposed 
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pressure  field,  the  usual  steady  state  boundary  layer  equations  are  unstable 
when  solved  as  an  initial  value  problem  in  space  in  regions  of  reversed  flow. 
However,  the  equations  can  be  marched  in  space  by  suppressing  the  streamwise 
convection  term  in  the  separated  region  (Ref.  6).  Although  this  approximation 
allows  the  solution  to  be  marched  through  separation,  the  approximation  becomes 
progressively  more  inaccurate  as  the  extent  of  the  separation  zone  or  the 
magnitude  of  either  the  normal  or  the  back  flow  velocities  become  large.  Thus, 
calculated  flow  details  which  may  be  important  (such  as  heat  transfer  at 
reattachment)  may  have  significant  error  when  separation  is  present  and  such  an 
interactive  analysis  is  used.  Other  schemes  have  and  are  being  developed  which 
solve  the  interaction  problem  without  encountering  an  instability  by  either 
changing  the  problem  to  initial  value  in  time  or  iteration  space.  Nonetheless 
the  resulting  solutions  still  retain  the  approximations  of  the  boundary  layer 
equations  and  the  inviscid  flow.  However,  as  with  invlscid  flow  solutions, 
combined  viscous  and  inviscid  solutions  remain  a  valuable  tool  for  those 
classes  of  cascade  problems  where  the  approximations  adopted  are  valid. 

The  final  procedure  currently  available  is  the  solution  of  full 
ensemble-averaged  Navier-Stokes  equations.  Such  an  analysis  has  been  applied 
to  a  variety  of  cascade  flow  fields  by  Shamroth,  Gibeling  and  McDonald  (Ref.  7) 
Shamroth,  McDonald  and  Briley  (Refs.  5,  8-10)  and  by  Shamroth  and  McDonald 
(Ref.  4).  The  use  of  the  full  Navier-Stokes  equations  for  the  cascade  problem 


6.  Rehyner,  T.  and  Flugge-Lotz,  I.:  The  Interaction  of  Shock  Waves  with  a 
Laminar  Boundary  Layer.  International  Journal  of  Nonlinear  Mechanics, 

Vol .  3,  1968. 

7.  Shamroth,  S.J.,  Gibeling,  H.J.  and  McDonald,  H. :  A  Navier-Stokes  Solution 
of  Laminar  and  Turbulent  Flow  Through  a  Cascade  of  Airfoils.  AIAA  Paper 
No.  80-1426,  1980.  (Also,  SRA  Report  R79-920004-F,  1982.) 

8.  Shamroth,  S.J.,  McDonald,  H.  and  Briley,  W.R. :  A  Navier-Stokes  Solution 
for  Transonic  Flow  Through  a  Cascade.  SRA  Report  R82-920007-F,  1982. 

9.  Shamroth,  S.J.,  McDonald,  H.  and  Briley,  W.R.:  Application  of  a  Navier- 
Stokes  Analysis  to  Transonic  Cascade  Flow  Fields.  ASME  Paper  82-GT-235, 
1982. 

10.  McDonald,  H.,  Shamroth,  S.J.  and  Briley,  W.R.:  Transonic  Flows  with 

Viscous  Effects,  Transonic  Shock  and  Multi -Dimensional  Flows:  Advances  in 
Scientific  Computing.  Academic  Press,  New  York,  1982. 
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allows  use  of  a  single  set  of  equations  for  the  entire  flow  field  and  thus 
removes  the  need  for  an  interaction  analysis  to  couple  different  equation 
descriptions  for  different  flow  regions.  The  analysis  simultaneously  predicts 
both  the  blade  pressure  distribution  and  viscous  and  heat  transfer  effects. 

The  initial  effort  in  the  present  ongoing  program  is  detailed  in  Ref.  7 
where  a  constructive  cascade  coordinate  system  was  combined  with  a 
Navier-Stok.es  calculation  procedure  to  compute  subsonic  flow  fields  in  a  simple 
cascade  of  unstaggered  NACA  0012  airfoils.  The  coordinate  system  generation 
process  was  based  upon  that  of  Eiseman  (Ref.  11).  The  calculation  procedure 
used  was  the  linearized  block  implicit  (LBI)  method  of  Briley  and  McDonald 
(Ref.  12)  which  obtains  a  solution  by  marching  the  assumed  initial  flow  field 
in  time  until  a  steady  state  is  reached.  Although  the  cases  chosen  were 
geometrically  simple  cascade  flows,  the  convergence  of  the  solution  to  steady 
state  in  a  relatively  few  number  of  time  steps  (~150)  showed  the  practicality 
of  the  approach. 

The  approach  was  extended  to  more  realistic  cascade  flow  fields  in  Ref.  8 
where  calculations  for  a  cascade  of  Sanz  airfoils  were  made.  The  coordinate 
system  generation  code  developed  by  Eiseman  was  modified  by  Kim  and  Shamroth 
(Ref.  13)  to  allow  specification  of  more  general  cascades  than  were  previously 
allowed.  This  new  generalized  code  was  used  to  generate  the  Sanz  cascade 
coordinate  system.  In  addition,  under  this  effort,  a  study  was  made  of 
possible  artificial  dissipation  formulations.  Based  upon  results  obtained  for 
a  model  one-dimensional  flow  problem  containing  a  shock  wave,  a  second  order 
artificial  dissipation  approach  was  judged  to  be  an  effective  method  of 
suppressing  spatial  oscillations  resulting  from  a  central  difference 
representation  of  derivatives.  This  same  technique  was  then  applied  to 
transonic  flow  through  a  cascade  of  Sanz  airfoils.  It  was  demonstrated  that 


11.  Eiseman,  P.R. :  A  Coordinate  System  for  a  Viscous  Transonic  Cascade 
Analysis.  Journal  of  Computational  Physics,  Vol.  ^6,  March  1978, 
pp.  307-338. 

12.  Briley,  W.R.  and  McDonald,  H. :  Solution  of  the  Multi-Dimensional 
Compressible  Navier-Stokes  Equations  by  a  Generalized  Implicit  Method. 
Journal  of  Computational  Physics,  Vol.  24,  1977. 

13.  Kim,  Y.— N.  and  Shamroth,  S.J.:  Revised  Coordinate  Generation  Program  for  a 
Cascade  of  Arbitrary  Shaped  Airfoils.  SRA  Report  81-1,  1981. 
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second  order  artificial  dissipation  could  be  added  to  the  governing  equations 
in  a  manner  which  would  suppress  spurious  spatial  oscillations  but  not 
significantly  contaminate  the  solution. 

Finally,  in  Ref.  4  the  same  procedure  was  used  to  calculate  flow  through 
turbine  and  compressor  cascade  where  experimental  data  was  available.  The 
turbine  cascade  chosen  was  that  of  Turner  (Ref.  14)  and  the  compressor  cascade 
that  of  Stephens  and  Hobbs  (Refs.  15  and  16).  Both  subsonic  and  transonic  flow 
conditions  were  considered  and  comparisons  between  calculation  and  measurement 
were  made  for  both  pressure  distributions  and  velocity  profiles.  In  general 
the  comparison  for  blade  pressure  distribution  was  judged  to  be  good  (see  Ref. 

4  for  details),  however  some  discrepancy  appeared  in  the  boundary  layer 
profiles.  The  present  effort  focuses  upon  a  further  study  of  the  boundary 
layer  profile  comparison  as  well  as  consideration  of  boundary  conditions  to  be 
applied  for  blade  movement  corresponding  to  small  forced  vibrations. 


14.  Turner,  A.B.:  Local  Heat  Transfer  Measurements  on  a  Gas  Turbine  Blade. 
Journal  of  Mechanical  Engineering  Sciences,  Vol.  13,  1971. 

15.  Stephens,  H.E.  and  Hobbs,  D.E.:  Design  and  Performance  of  Supercritical 
Airfoils  for  Axial  Flow  Compressors.  Pratt  and  Whitney  Aircraft 
Report  FR11455,  1979. 

16.  Hobbs,  D.E.,  Wagner,  J.H.,  Dannenhoffer ,  J.F.,  Dring,  R.P.:  Wake 
Experiment  and  Modelling  for  Fore  and  Aft— loaded  Compressor  Cascade. 
Pratt  and  Whitney  Aircraft  Report  FR13514,  1980. 
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II.  ANALYSIS 


The  present  analysis  is  based  upon  a  solution  of  the  ensemble-averaged 
Navier-Stokes  equations  using  the  linearized  block  implicit  method  of  Briley 
and  McDonald  (Ref.  12).  The  equations  are  solved  in  a  constructive  coordinate 
system  (Ref.  4)  with  density  and  the  Cartesian  velocity  components  being  taken 
as  dependent  variables.  The  application  of  the  LBI  method  to  the  cascade  flow 
field  problem  has  been  discussed  in  some  detail  in  Refs.  4,  and  7-10.  However, 
for  completeness  it  will  be  repeated  here  along  with  a  brief  discussion  of  the 
coordinate  system  and  governing  equations. 

Coordinate  System 

An  important  component  of  the  cascade  analysis  is  the  cascade  coordinate 
system.  Any  coordinate  system  used  in  the  analysis  should  satisfy  conditions 
of  (i)  generality,  (ii)  smoothness,  (iii)  resolvability  and  (iv)  allow  easy 
application  of  boundary  conditions.  Obviously,  a  coordinate  system  must  be 
sufficiently  general  to  allow  application  to  a  wide  class  of  problems  of 
interest  if  it  is  to  be  practical.  The  metric  data  associated  with  the  coordi¬ 
nate  system  must  be  sufficiently  smooth  so  that  the  variation  from  grid  point 
to  grid  point  does  not  lead  to  a  numerical  solution  dominated  by  metric  coef¬ 
ficient  truncation  error.  It  should  be  noted  that  this  requirement  differs 
from  the  requirement  of  the  existence  of  a  specified  number  of  transformation 
derivatives.  The  coordinate  system  must  resolve  flow  regions  where  rapid  flow 
field  changes  occur.  Finally,  coordinates  should  allow  accurate  implementation 
of  boundary  conditions;  for  the  cascade  this  requirement  is  equivalent  to  the 
requirement  that  the  metric  coefficients  be  continuous  across  the  periodic 
lines  where  periodic  boundary  conditions  are  to  be  applied. 

To  date,  several  types  of  coordinate  systems  are  available.  These  include 
(i)  solutions  based  upon  a  conformal  transformation,  (ii)  solutions  based  upon 
solution  of  a  Poisson  equation  (e.g.  Ref.  17),  and  constructive  systems.  The 
present  approach  uses  a  constructive  system  based  originally  upon  the  approach 
of  Eiseman  (Ref.  11)  which  has  been  revised  by  Kim  and  Shamroth  (Ref.  13).  A 

17.  Thompson,  J.F.,  Thames,  F.C.  and  Mastin,  C.W. :  Boundary  Fitted  Curvi¬ 
linear  Coordinate  Systems  for  Solution  of  Partial  Differential  Equations 
on  Fields  Containing  Any  Number  of  Arbitrary  Two-Dimensional  Bodies. 

NASA  CR-2729,  July  1977. 


6 


computer  generated  plot  of  the  coordinate  system  for  the  Jose  Sanz  diffusion 
cascade  is  shown  in  Fig.  1.  As  can  be  seen  in  Fig.  1,  the  coordinate  systems 
consist  of  two  sets  of  curves;  the  C  =  constant  curves  such  as  FG  or  HI  and  the 
n  =  constant  curves  such  as  ABCD  or  A'ED'.  In  constructing  the  coordinate 
system  care  must  be  taken  that  metric  data  remains  smooth  from  grid  point  to 
grid  point,  and  adequate  resolution  is  obtained  both  near  the  blade  surface  and 
in  the  leading  edge  region.  Details  of  the  construction  process  are  given  in 
Appendix  A. 


Governing  Equations 

The  present  effort  solves  the  time-dependent  compressible  Navier— Stokes 
equations  to  predict  the  cascade  flow  field.  If  the  computational  spatial 
coordinates  are  £  and  p  where 

€m  t)  *  ^(x  ,y,  t)  r«t  (D 

then  the  continuity  equation,  the  x-component  of  the  momentum  equation  and  the 
y-component  of  the  momentum  equation  are  written  as 


dW  ^  dW  *  dF  ^  dG  dW  dF  dG 

'TT  +  (<  W  +  W  *  IT  +  *><Tj  +  + 


d  -  fxljy  -  (yVx 


In  Eqs .  (1-3)  x  and  y  are  Cartesian  coordinates,  t  is  time,  u  and  v  are 
velocity  components,  p  is  density,  p  is  pressure,  and  x-^j  is  the  stress 
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tensor  and  Re  is  the  Reynolds  number.  This  formulation  is  termed  the 
quasi-linear  form  and  has  been  used  successfully  for  a  number  of  cascade  and 
airfoil  calculations  [7-10,  18,  19]  for  both  laminar  and  turbulent  flow  in  the 
subsonic  and  transonic  regimes. 

The  dependent  variables  chosen  for  the  present  formulation  are  the 
density,  p,  and  the  velocity  components,  u  and  w.  Although  the  code  does 
contain  an  energy  equation  and  calculations  have  been  made  with  an  energy 
equation,  most  calculations  have  been  run  with  the  assumption  T°,  the 
stagnation  temperature,  equals  a  constant.  With  this  assumption,  the  pressure 
is  related  to  the  velocity  and  density  by 


P 


pR 


\ 

2Cp  ) 


(4) 


Numerical  Procedure 


The  numerical  procedure  used  to  solve  the  governing  equations  is  a 
consistently  split  linearized  block  implicit  (LBI)  scheme  originally  developed 
by  Briley  and  McDonald  (Ref.  12).  A  conceptually  similar  scheme  has  been 
developed  for  two-dimensional  MHD  problems  by  Lindemuth  and  Killeen  (Ref.  20). 
The  procedure  is  discussed  in  detail  in  Refs.  12  and  21.  The  method  can  be 
briefly  outlined  as  follows :  the  governing  equations  are  replaced  by  an 
implicit  time  difference  approximation,  optionally  a  backward  difference  or 


18.  Shamroth,  S.J.  and  Gibeling,  H.J.:  A  Compressible  Solution  of  the  Navier- 
Stokes  Equations  for  Turbulent  Flow  about  an  Airfoil.  NASA  CR-3183,  1979. 
(Also,  AIAA  Paper  79-1543). 

19.  Shamroth,  S.J.  and  Gibeling,  H.J.:  Analysis  of  Turbulent  Flow  about  an 
Isolated  Airfoil  Using  a  Time-Dependent  Navier-Stokes  Procedure.  AGARD  CP 
296,  1980. 

20.  Lindemuth,  I.  and  Killeen,  J.:  Alternating  Direction  Implilcit  Techniques 
for  Two-Dimensional  Magnetohydrodynamic  Calculations.  Journal  of 
Computational  Physics,  13,  1973. 

21.  Briley,  W.R.,  and  McDonald,  H.:  On  the  Structure  and  Use  of  Linearized 
Block  Implicit  Schemes.  Journal  of  Computational  Physics,  Vol.  34,  1980. 
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Crank-Nicolson  scheme.  Terms  involving  nonlinearities  at  the  implicit  time 
level  are  linearized  by  Taylor  expansion  in  time  about  the  solution  at  the 
known  time  level,  and  spatial  difference  approximations  are  introduced.  The 
result  is  a  system  of  multi -dimensional  coupled  (but  linear)  difference 
equations  for  the  dependent  variables  at  the  unknown  or  implicit  time  level . 

To  solve  these  difference  equations,  the  Douglas-Gunn  (Ref.  22)  procedure  for 
generating  alternating-direction  implicit  (ADI)  schemes  as  perturbations  of 
fundamental  implicit  difference  schemes  is  introduced  in  its  natural  extension 
to  systems  of  partial  differential  equations.  This  technique  leads  to  systems 
of  coupled  linear  difference  equations  having  narrow  block-banded  matrix 
structures  which  can  be  solved  efficiently  by  standard  block-elimination 
methods. 

The  method  centers  around  the  use  of  a  formal  linearization  technique 
adapted  for  the  integration  of  initial -value  problems.  The  linearization 
technique,  which  requires  an  implicit  solution  procedure,  permits  the  solution 
of  coupled  nonlinear  equations  in  one  space  dimension  (to  the  requisite  degree 
of  accuracy)  by  a  one-step  noniterative  scheme.  Since  no  iteration  is  required 
to  compute  the  solution  for  a  single  time  step,  and  since  only  moderate  effort 
is  required  for  solution  of  the  implicit  difference  equations,  the  method  is 
computationally  efficient;  this  efficiency  is  retained  for  multi-dimensional 
problems  by  using  what  might  be  termed  block  ADI  techniques.  The  method  is 
also  economical  in  terms  of  computer  storage,  in  its  present  form  requiring 
only  two  time-levels  of  storage  for  each  dependent  variable.  Furthermore,  the 
block  ADI  technique  reduces  multi-dimensional  problems  to  sequences  of 
calculations  which  are  one— dimensional  in  the  sense  that  easily— solved  narrow 
block-banded  matrices  associated  with  one— dimensional  rows  of  grid  points  are 
produced.  A  more  detailed  discussion  of  the  solution  procedure  as  discussed  by 
Briley,  Buggeln  and  McDonald  (Ref.  23)  is  given  in  the  Appendix  B. 

Artificial  Dissipation 

Since  the  calculations  of  interest  are  often  at  high  Reynolds  numbers 

22.  Douglas,  J.  and  Gunn,  J.E.:  A  General  Formulation  of  Alternating 
Direction  Methods.  Numerische  Math.,  Vol.  6,  1964,  pp.  428-453. 

23.  Briley,  W.R.,  Buggeln,  R.C.  and  McDonald,  H.:  Computation  of  Laminar 
and  Turbulent  Flow  in  90  Degree  Square  Duct  and  Pipe  Bends  Using  the 
Navier-Stokes  Equations.  SRA  Report  R82-920009-F,  1982. 
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typical  of  normal  turbomachinery  applications,  it  is  necessary  to  add 
"artificial  dissipation"  terms  to  suppress  spatial  oscillations  associated  with 
central  spatial  differences  approximations.  This  can  be  done  via  a  dissipative 
spatial  difference  formulation  (e.g.,  one-sided  difference  approximations  for 
first  derivatives)  or  by  explicitly  adding  an  additional  dissipative  type 
term.  For  the  Navier-Stokes  equations,  the  present  authors  favor  the  latter 
approach  since  when  an  additional  term  is  explicitly  added,  the  physical 
approximation  being  made  is  usually  clearer  than  when  dissipative  mechanisms 
are  contained  within  numerical  truncation  errors,  and  further,  explicit 
addition  of  an  artificial  dissipation  term  allows  greater  control  over  the 
amount  of  non-physical  dissipation  being  added.  Obviously,  the  most  desirable 
technique  would  add  only  enough  dissipative  mechanism  to  suppress  oscillations 
without  deteriorating  solution  accuracy.  Various  methods  of  adding  artificial 
dissipation  were  investigated  in  Ref.  8,  and  these  were  evaluated  in  the 
context  of  a  model  one-dimensional  problem  containing  a  shock  with  a  known 
analytic  solution  (one-dimensional  flow  with  heat  transfer).  The  methods  which 
were  considered  included  second-order  dissipation,  fourth-order  dissipation  and 
pressure  dissipation  techniques. 

As  a  result  of  this  investigation,  it  was  concluded  that  a  second-order 
anisotropic  artificial  dissipation  formulation  suppressed  spatial  oscillations 
without  impacting  adversely  on  accuracy  and  could  be  used  to  capture  shocks 
successfully.  In  this  formulation,  the  terms 


P 


n-i 


ay2 


are  added  to  the  governing  equations  where  <f>  =  u,  v  and  p  for  the  x -momentum, 
y-momentum  and  continuity  equations,  respectively.  The  exponent  n  is  zero  for 
the  continuity  equation  and  unity  for  the  momentum  equations.  The  dissipation 
coefficient  dx  is  determined  as  follows.  The  general  equation  has  an 
x-direction  convective  term  of  the  form  a3<J>/3x  and  an  x-direction  diffusion 
term  of  the  form  3(b3<J>/3x)/3x.  The  diffusive  term  is  expanded 
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d(b  dcfi/ dx)/dx  =  b  d2<£/dx2  +  db/dx  dcfc/dx 


(5) 


and  then  a  local  cell  Reynolds  number  Re^x  is  defined  for  the  x-direction  by 


ReAx  =  |a  ~  db/3x|Ax/b 


(6) 


where  b  is  the  total  or  effective  viscosity  including  both  laminar  and 
turbulent  contributions,  and  Ax  is  the  grid  spacing.  The  dissipation  coef¬ 
ficient  dx  is  non-negative  and  is  chosen  as  the  larger  of  zero  and  the  local 
quantity  b  (<JxRe^x-l).  The  dissipation  parameter  is  a  specified 
constant  and  represents  the  inverse  of  the  cell  Reynolds  number  below  which  no 
artificial  dissipation  is  added.  The  dissipation  coefficient  dy  is  evaluated 
in  an  analogous  manner  and  is  based  on  the  local  cell  Reynolds  number  Re^y 
and  grid  spacing  Ay  for  the  y-direction  and  the  specified  parameter  Oy.  It 
should  be  noted  that  recently  calculations  have  been  run  with  artificial  dissi¬ 
pation  added  in  the  conservative  form  3(  pn“ld  3$/  3x)/ 3x  and  no 
significant  difference  between  the  forms  was  noted  for  the  particular  flows 
examined. 

The  question  arises  as  to  the  values  of  ax  and  Oy  which  should  be 
chosen.  This  was  assessed  both  through  the  model  problem  (Ref.  8),  and  through 
calculations  for  a  Jose  Sanz  compressor  cascade  (Refs.  4  and  8)*  These  results 
indicated  that  values  of  a  =  .5  which  correspond  to  a  cell  Reynolds  number  2 
limitation  would  severely  damp  physical  variations.  However,  when  a  was  set 
in  the  range  .025  -  cr  ^  0.05,  which  correspond  to  a  cell  Reynolds  number  range 
between  40  and  20,  spurious  spatial  oscillations  were  damped  with  no 
significant  change  in  the  calculated  results  as  a  was  varied  in  this  range. 
Further,  as  discussed  in  Refs.  4  and  7-10,  the  results  obtained  showed  good 
agreement  with  data.  This  has  since  been  confirmed  at  several  other  studies  at 
Scientific  Research  Associates  such  as  two-  and  three-dimensional  transonic 
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nozzle  flows  (Ref,  24)  where  a  maximum  acceptable  value  of  a  =  0,10  has  been 
noted  for  most  problems.  Therefore,  based  upon  this  past  experience, 
second-order  damping  is  applied  with  a  taken  as  0.05. 

B  ound  ar  y  Cond i t i ons 

The  authors 1  experience  in  solving  Navier-Stokes  equations  has  indicated 
the  important  role  boundary  conditions  play  in  determining  accurate  solutions 
and  rapid  numerical  convergence.  The  boundary  conditions  used  in  the  present 
calculations  follow  the  suggestion  of  Briley  and  McDonald  (Ref.  25)  which 
specifies  upstream  total  pressure  and  downstream  static  pressure  conditions. 
For  the  cascade  system  shown  in  Fig.  1,  AB  and  CD  are  periodic  boundaries  and 
periodic  conditions  are  set  here. 

Specification  of  upstream  and  downstream  conditions  is  somewhat  more 
difficult.  For  an  isolated  cascade,  boundary  conditions  for  the  differential 
equations  may  be  known  at  both  upstream  infinity  and  downstream  infinity. 
However,  since  computation  economics  argues  for  placing  grid  points  in  the 
vicinity  of  the  cascade  and  minimizing  the  number  of  grid  points  far  from  the 
cascade,  the  upstream  and  downstream  computational  boundaries  should  be  set  as 
close  to  the  cascade  as  is  practical.  In  addition,  with  the  particular 
body-fitted  coordinates  used,  as  the  upstream  boundary  moves  further  upstream, 
the  angle  between  pseudo-radial  and  pseudo-azimuthal  coordinate  lines  becomes 
smaller.  Decreasing  the  coordinate  angle  causes  the  coordinate  system  to 


24.  Liu,  N.S.,  Shamroth,  S.J.  and  McDonald,  H.:  Numerical  Solution  of  the 
Navier-Stokes  Equations  for  Compressible  Turbulent  Two/Three  Dimensional 
Flows  in  the  Terminal  Shock  Region  of  an  Inlet/Diffuser.  AIAA  Paper 
83-1892,  1983. 

25.  Briley,  W.R.  and  McDonald,  H. :  Computation  of  Three-Dimensional  Horseshoe 
Vortex  Flow  Using  the  Navier-Stokes  Equations.  Seventh  International 
Conference  on  Numerical  Methods  in  Fluid  Dynamics,  1980. 


12 


to  become  less  well-conditioned,  increases  truncation  error  (Ref,  26),  and 
increases  the  role  of  cross-derivative  terms  in  the  equations.  All  of  these 
characteristics  could  be  detrimental  to  the  present  numerical  procedure  and, 
therefore,  they  also  argue  for  placing  the  upstream  boundary  as  close  to  the 
cascade  as  possible.  However,  when  the  upstream  boundary  is  placed  close  to 
the  cascade,  most  flow  function  conditions  on  the  boundary  will  not  be  known, 
since  these  will  have  been  changed  from  values  at  infinity  by  the  presence  of 
the  cascade. 

In  the  present  approach,  the  suggestion  of  Ref.  25  is  followed  which  sets 
total  pressure  on  boundary  BC  (see  Fig.  1).  Unless  boundary  BC  is  very  far 
upstream,  the  flow  velocity  along  this  boundary  will  not  be  equal  to  the 
velocity  at  upstream  infinity  since  some  inviscid  deceleration  will  have 
occurred.  However,  as  long  as  the  boundary  is  upstream  of  the  region  of  any 
significant  viscous  or  shock  phenomena,  the  total  pressure  on  this  boundary 
will  be  equal  to  the  total  pressure  at  upstream  infinity.  Hence,  total 
pressure  is  an  appropriate  boundary  condition  realistically  modeling  the 
desired  flow  condition.  In  addition  to  specifying  upstream  total  pressure,  it 
is  necessary  to  specify  the  inlet  flow  angle.  In  the  present  calculation,  a 
value  was  assumed  constant  on  the  upstream  boundary  at  a  specified  value.  The 
third  condition  set  on  the  upstream  boundary  concerns  the  density  and  a  zero 
density  derivative  at  this  boundary  was  specified  as  a  numerical  treatment  of 
the  boundary.  The  downstream  boundary  was  treated  by  setting  a  constant  static 
pressure  as  a  boundary  condition,  and  by  setting  second  derivatives  of  both 
velocity  components  equal  to  zero  at  this  location.  In  the  present 
application,  a  constant  static  pressure  was  set  at  downstream  infinity,  and 
hence  it  is  assumed  that  the  downstream  boundary  is  located  in  a  region  where 
pressure  is  uniform. 

Both  the  upstrem  and  downstream  boundaries  have  boundary  conditions 
associated  with  them  which  are  nonlinear  functions  of  the  dependent  variables. 
These  are  the  specifications  of  total  pressure  on  the  upstream  boundary  and 
static  pressure  on  the  downstream  boundary.  These  nonlinear  boudnary 
conditions  are  linearized  in  the  same  manner  as  the  governing  equations  (via  a 
Taylor  expansion  of  the  dependent  variable  in  time),  and  then  solved  implicitly 

26.  Mastin,  C.W.:  Error  Induced  by  Coordinate  Systems,  Numerical  Grid 

Generation,  J.F.  Thompson,  Ed.,  Elsevier  Publishing, _New  York,  1982. 
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along  with  the  interior  point  equations.  The  final  boundary  conditions  to  be 
considered  are  the  conditions  along  the  blade  surface.  Here  no-slip  and  no 
through-flow  conditions  were  applied  leading  to  a  specification  of  zero 
velocity  on  the  surface.  An  inviscid  transverse  momentum  equation  was  applied 
on  the  surface  leading  to  a  boundary  condition  approximation  of  zero  transverse 
pressure  gradient  being  applied. 

III.  CURRENT  EFFORTS 

Current  efforts  pursued  under  the  present  contract  consisted  of  three 
parts:  computer  code  speed-up,  turbulence  model  studies  and  consideration  of 
time-dependent,  forced  vibration-type  boundary  conditions.  These  are  now 
discussed  in  detail. 

Code  Efficiency 

The  original  cascade  computer  code  was  a  very  general  code  which  was 
written  for  maximum  flexibility.  In  this  regard,  equations,  boundary 
conditions,  dependent  variables,  coordinate  systems,  difference  schemes,  etc. 
could  be  easily  changed  and  the  choice  of  equations  and  boundary  conditions 
were  made  by  the  user.  For  example,  user  specified  input  determined  if  an 
energy  equation  was  solved  or  constant  total  enthalpy  assumed,  the  choice  of 
two-  or  three-dimensional  flow  was  made  via  input  and  the  choice  of 
constructive  versus  quasi-linear  equations  form  was  made  by  input.  In 
addition,  several  options  such  as  reacting  flow,  two-phase  flow,  cylindrical 
polar  geometry,  etc.  which  are  not  relevant  to  the  cascade  problem  were 
available.  Although  this  generality  is  very  desirable  from  the  viewpoint  of 
deck  development  and  deck  flexibility,  it  does  require  a  price  in  terms  of 
computer  run  time.  Therefore,  as  a  first  item  under  the  present  effort  the 
existing  code  was  modified  to  increase  deck  efficiency. 

The  modifications  still  allow  considerable  generality  such  as  choice  of 
equations,  choice  of  boundary  conditions,  arbtrary  geometry,  conservative  or 
quasi-linear  equation  form,  etc.  However,  certain  portions  of  the  code  not 
required  for  cascade  calculations  such  as  cylindrical  polar  coordinates  and 
reacting  flow  were  eliminated.  In  addition,  the  original  matrix  inverter  was 
replaced  by  a  more  efficient  procedure  and  the  computational  coordinates  were 
assumed  to  be  equally  spaced.  It  should  be  noted  that  since  the  computational 


14 


coordinates  are  general  non-physical  coordinates,  see  Eq.  (1),  this  does  not 
place  a  restriction  on  nonuniform  meshes  in  physical  space  and,  in  fact,  all 
cases  run  to  date  were  run  with  equally  spaced  computational  coordinates. 

The  resulting  modification  decreased  computer  run  time  by  approximately  55 
per  cent  and,  therefore,  allows  considerable  additional  computer  runs  for  a 
given  computer  budget.  The  new  code  requires  approximately  0.0034  secs  per 
grid  point  per  time  step  for  a  code  operating  with  an  out-of-core  option.  If 
the  code  were  run  completely  in  core,  this  time  would  reduce  to  approximately 
0.0025  secs.  It  is  estimated  that  an  additional  factor  of  2  saving  could  be 
obtained  without  loss  of  generality  with  further  code  restructuring. 

Obviously,  even  further  savings  are  available  with  code  vector ization  and/or 
modifications  restricting  generality. 


Turbulence  Model  Modifications 


One  major  focus  of  the  present  effort  was  further  consideration,  and 
possible  modification,  of  turbulence  models  in  the  existing  computer  program. 

At  the  beginning  of  the  present  phase,  the  code  allowed  for  either  of  two 
turbulence  models;  these  were  a  mixing  length  model  and  a  turbulence  energy  - 
algebraic  length  scale  model. 

In  the  mixing  length  model,  the  turbulent  viscosity  is  related  to  the  mean 
strain  via  a  mixing  length,  £,  such  that 


1/2 


(7) 


where  pT  is  the  turbulent  viscosity,  p  is  the  density,  i  is  the  mixing 
length,  u^  Is  the  velocity  component  and  x-£  is  the  Cartesian 


direction.  Summation  is  implied  for  the  repeated  indices.  The  question  now 
arises  as  to  specification  of  £.  For  the  region  upstream  of  the  trailing  edge 
the  mixing  length  is  specified  in  the  usual  boundary  layer  manner;  i.e. 


I 


i  <  l 


(8) 
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where  k  is  the  von  Karman  constant  and  y+  is  the  dimensionless  normal 
coordinate,  yuT/v.  In  boundary  layer  analysis  Vax  *s  usually  taken  as 
0.096  where  6  is  the  boundary  layer  thickness  taken  as  the  location  where 
u/ue  =  0.99.  However,  this  definition  of  6  assumes  the  existence  of  an  outer 
flow  where  the  velocity  ue  is  independent  of  distance  from  the  wall  at  a 
given  streamwise  station,  i.e.,  it  assumes  ue  is  only  a  function  of  the 
streamwise  coordinate.  Although  a  boundary  layer  calculation  will  yield 
solutions  in  which  u  approaches  ue  asymptotically  at  distances  far  from  the 
solid  no-slip  surface,  Navier-Stokes  solutions  for  cascade  flow  fields  do  not 
in  general  predict  a  region  where  u  asymptotes  to  a  constant  value. 
Furthermore,  measurements  of  the  flow  also  show  no  such  region  to  exist  in 
general.  Obviously,  a  proper  choice  of  6  for  the  Navier-Stokes  cascade 
analysis  is  not  straight  forward.  Calculations  made  in  Refs.  4  and  7-9  have 
set  the  boundary  layer  thickness  by  first  determining  u^^.,  the  maximum 
streamwise  velocity,  at  a  given  station  and  then  setting  6  via 


s  = 


2.0  y 

(u/umax=  k|> 


(9) 


;i.e.,  6  was  taken  as  twice  the  distance  for  which  u/umax  =  kj.  Two  values 
of  k^  were  used  in  the  previous  efforts;  these  were  0.80  and  0.95.  Although 
the  predicted  pressure  distribution  was  relatively  insensitive  to  choice,  the 
boundary  layer  development  did  depend  upon  the  choice  of  k^. 

The  model  used  in  the  wake  is  also  a  mixing  length  model  in  which  the 
mixing  length  was  made  proportional  to  the  wake  height,  6,  and  a  linear  of 
growth  of  6  with  distance  was  assumed  based  upon  the  classical  free  jet 
boundary  results  (e.g.  [27]).  With  the  free  jet  boundary  growth  assumption 


^  =  (SpS  +  Sss)  +  ( .2)(x  xTE) 


(10) 


where  5pS  and  6gs  are  the  pressure  and  suction  surface  trailing  edge 
boundary  layer  thickness  and  xte  is  the  trailing  edge  location.  The  mixing 
length,  £,  was  taken  as  0.26. 


27.  Schlichting,  H. :  Boundary  Layer  Theory,  McGraw  Hill,  New  York, 
1960. 
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Using  this  formulation,  calculated  profiles  were  compared  with  measured 
profiles  for  the  data  of  Hobbs,  Wagner,  Dannenhoffer  and  Dring  (Ref.  16).  In 
Ref.  16,  boundary  layer  profiles  were  measured  on  both  the  suction  and  pressure 
surfaces  in  the  vicinity  of  the  blade  trailing  edge  region.  In  this  regard,  it 
should  be  noted  that  the  resulting  boundary  layers  on  both  surfaces  are 
subjected  to  strong  favorable  pressure  gradients  followed  by  adverse  pressure 
gradients  with  these  effects  being  considerably  more  pronounced  on  the  suction 
surface.  Therefore,  the  boundary  layer  at  the  trailing  edge  has  a  history  of 
both  strong  favorable  and  adverse  pressure  gradients  and  its  profile  depends 
both  upon  this  history  and  the  transition  point  location. 

The  comparisons  of  Ref.  4  showed  very  good  agreement  between  the  measured 
profiles  and  the  calculated  profiles  on  the^pressure  surface.  On  the  suction 
surface  agreement,  although  acceptable,  was  not  as  good.  In  addition,  the 
choice  of  the  variable,  kj ,  (see  Eq.  (9))  could  significantly  change  the 
calculated  profile.  Therefore,  a  more  detailed  study  of  the  turbulence  model 
was  warranted. 

As  a  first  step  in  the  re-evaluation,  the  code  was  used  to  calculate  a 
constant  pressure  turbulent  boundary  and  the  resulting  profile  is  compared  to 
data  compilations  in  Fig.  2.  As  can  be  seen,  the  predicted  profile  is  in  good 
agreement  with  the  compiled  data.  In  this  calculation 

8  =  2.0y 

(u/umax=0-8)  <U> 

These  results  in  conjunction  with  the  boundary  layer  profiles  presented  in 
Ref.  4  indicate  reasonable  agreement  between  calculation  and  experiment  can  be 
obtained  using  the  mixing  length  model  with  a  proper  choice  of  6.  Obviously, 
these  results  are  preliminary  and  not  conclusive.  To  be  made  conclusive, 
considerably  more  comparisons  with  data  would  be  required. 

The  major  obstacle  in  using  the  mixing  length  model  resides  in  the 
specification  of  the  length  scale  6,  and  it  may  be  advantageous  to  utilize  a 
turbulence  model  which  does  not  require  specification  for  a  length  scale.  One 
model  which  fits  this  requirement  is  the  so-called  k-e  model  which  solves 
governing  equations  for  the  turbulence  energy,  k,  and  the  turbulent 
dissipation,  e,  and  then  calculates  yT  from  these  quantities.  The  equations 
governing  the  development  of  k  and  e  have  been  given  by  several  authors 
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(e.g.  Launder  and  Spalding  Ref.  28).  One  form  in  which  the  equations  may  be 
expressed  is 
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Mt  =  Pcp*Z /€ 


(14) 


where  D,  E,  Cj,  C2,  Cy ,  and  ae  are  functions  which  vary  from 

particuar  model  to  model.  A  recent  comparison  of  several  k-e  models  has  been 

given  by  Patel,  Rodi  and  Scheuerer  (Ref.  29).  In  the  present  analysis  the 


28.  Launder,  B.E.  and  Spalding,  D.B.:  The  Numerical  Computation  of  Turbulent 
Flows,  Computer  Methods  in  Applied  Mechanics  and  Engineering,  Vol.  3, 

1974,  pp.  269-287. 

29.  Patel,  V.C.,  Rodi,  W.  and  Scheurer,  G.:  Evaluation  of  Turbulence  Models 
for  Near  Wall  and  Low  Reynolds  Number  Flows.  Third  Symposium  on  Turbulent 
Shear  Flow.  University  of  California,  Davis,  1981. 
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following  functional  forms  were  set 


cr€  b  1.3 


o-k  =  1.0 
C,  =  1.43 

=  0.09  exp  [-  2.5/( I  +  Rt/50)] 
C2  =  1.92  [l.O  -  0.3  exp  (-R|)] 


D 


E 


■2>j. 


akl/2  akl/2 

aXj  ax. 


2^t 


a2 


U, 


d*k  \ 


(15) 


This  corresponds  to  equations  given  by  Launder  and  Spalding  in  Ref.  28. 

Further  discussion  of  the  equations  is  given  in  some  detail  in  Refs.  28  and  29. 

A  similar  Navier-Stokes  code  using  the  same  equations  and  the  same 
algorithm  has  been  run  at  SRA  under  an  Army  Research  Office  program  to 
calculate  a  transonic  shock  boundary  layer  interaction  in  a  constant  radius 
tube  (Ref.  30),  and  these  results  are  repeated  here  since  they  demonstrate  a 


30.  McDonald,  H. ,  Roscoe,  D.V.,  Gibeling,  H.J.  and  Shamroth,  S.J.:  Progress 
Report  for  Contract  DAAG29-80-C-0082,  December  1981. 
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transonic  calculation.  The  case  considered  had  supersonic  inflow  in  a  straight 
tube  with  a  near  normal  shock  resulting  from  the  specified  back  pressure.  The 
calculation  was  run  to  correspond  with  data  of  Mateer,  Brosh  and  Viegas 
for  both  a  mixing  length  and  a  k-e  model  (Ref.  31).  A  comparison  of  calculated 
and  measured  pressure  distributions  is  given  in  Fig.  3.  As  can  be  seen, 
agreement  is  very  good.  A  comparison  for  the  predicted  and  measured  streamwise 
velocity  profiles  is  presented  in  Fig.  4.  Both  calculations  show  good 
agreement  with  data  although  the  mixing  length  calculation  is  in  slightly 
better  agreement  with  data  than  is  the  k-e  calculation.  These  results  give 
encouragement  for  the  use  of  the  k-e  model  in  non-simple  viscous  flow  fields 
where  a  mixing  length  application  may  not  suffice  and  give  impetus  to  the 
inclusion  of  thse  models  in  the  cascade  deck.  However,  as  discussed  in  Ref. 

29,  k-e  models  are  not  always  satisfactory  even  in  simple  boundary  layers  so 
the  approach  is  taken  with  some  caution. 

Under  the  present  effort,  Eqs.  (12)  and  (13)  have  been  incorporated  into 
the  Navier— Stokes  cascade  code  and  the  subsonic  cascade  configuration  of 
Hobbs,  et  al  (Ref.  16)  has  been  recalculated.  A  plot  of  predicted  and  measured 
blade  pressure  distribution  is  shown  in  Fig.  5  which  compares  measured  data 
with  mixing  length  and  k-e  calculations.  As  can  be  seen,  both  calculated 
pressure  distributions  are  nearly  the  same  and  agree  well  with  the  measured 
data. 

Predictions  of  the  boundary  layer  profile  are  given  in  Figs.  6-9. 

Although  the  calculations  presented  in  Ref.  4  give  good  agreement  with  data, 
they  were  run  for  a  value  of  kj  =  0.95  in  Eq.  (9)  which  appears  to  be 
unrealistically  large.  Therefore,  the  mixing  length  calculations  were  rerun 
for  k  =  0.90  which  appears  to  be  a  more  reasonable  value.  Calculations  were 
also  run  for  fully  turbulent  and  forced  transition  situations.  It  should  be 
recalled  that  when  the  mixing  length  model  is  used,  the  flow  is  considered 
turbulent  throughout  with  the  turbulent  viscosity  over  the  airfoil  set  by  Eqs. 
(7)  through  (11).  Since  the  maximum  length  scale  is  determined  by  the  boundary 
layer  thickness,  the  turbulent  viscosity  is  small  in  the  leading 

31.  Mateer,  G.G.,  Brock,  A.  and  Viegas,  J.R.:  A  Normal  Shock— Wave  Turbulent 

Boundary  Layer  Interaction  at  Transonic  Speeds.  AIAA  Paper  76-161,  1976. 
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edge  region  and  the  flow  is  essentially  laminar  here.  However,  a  turbulent 
viscosity  is  calculated  throughout;  this  calculation  is  termed  the  fully 
turbulent  calculation.  Alternatively,  the  calculation  can  be  run  in  a  mode 
where  the  turbulent  viscosity  is  set  to  zero  upstream  of  a  specified  streamwise 
location.  This  is  termed  the  forced  transition  model.  Similarly  when  the  k-e 
model  is  utilized,  it  predicts  a  turbulent  viscosity  throughout  the  flow.  In 
certain  regions,  this  viscosity  may  be  small  and  in  these  regions  the  flow  will 
be  laminar.  However,  a  calculation  can  be  made  whereby  the  turbulent  viscosity 
is  set  to  zero  upstream  of  a  specified  location.  This  latter  mode  is  termed 
the  forced  transition  mode. 

Predictions  of  the  pressure  surface  boundary  layer  are  given  in  Figs.  6 
and  7.  Figs.  6  and  7  compare  fully  turbulent  and  free  transition  models  for 
the  mixing  length  and  k-e  models  respectively.  As  can  be  seen,  agreement  in 
all  cases  is  good.  In  all  cases  of  forced  transition,  laminar  flow  was  imposed 
for  s/c  <  0.18  which  is  the  location  where  transition  steps  were  placed  in  the 
experiment. 

Predictions  of  the  suction  surface  boundary  layer  is  given  in  Figs.  8  and 
9.  The  mixing  length  calculations  are  given  in  Fig.  8  where  as  shown  the 
forced  transition  calculation  is  in  much  better  agreement  with  the  data  than  is 
the  fully  turbulent  calculation.  This  calculation,  in  conjunction  with  other 
calculations,  indicate  that  the  calculated  suction  surface  boundary  layer 
profile  is  sensitive  both  to  the  chosen  value  of  k^  and  the  assumption  of 
free  or  forced  transition.  The  major  effect  of  requiring  laminar  flow  for 
s/c  <  .18  was  to  maintain  a  thinner  boundary  layer  at  the  beginning  of  adverse 
pressure  gradient  regions.  This  effect  was  then  magnified  as  the  boundary 
layer  was  subjected  to  the  strong  adverse  pressure  gradient  and  the  results  are 
shown  in  Fig.  8. 

Calculations  made  with  the  k-e  model  are  given  in  Fig.  9.  As  can  be  seen, 
differences  again  appear  between  the  fully  turbulent  and  forced  transition 
model.  However,  the  effect  of  requiring  laminar  flow  for  s/c  <  .18,  is  not  as 
pronounced  in  this  case  and  both  calculations  predict  boundary  layers  which  are 
thicker  than  those  measured  experimentally.  It  should  be  noted  that  the  k-e 
equations  do  contain  low  Reynolds  number  effects.  These  are  incorporated  in 
the  model  via  the  functional  forms  Cy  and  C2  and  may  account  for  near  wall 
and/or  transitional  effects.  Although  these  models  have  been  used  with  some 
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success  for  relaminarizing  boundary  layers  (Ref.  32),  their  application  to 
boundary  layers  in  arbitrary  pressure  gradients  has  found  only  modest  success 
(e.g.  Ref.  29),  and  they  have  yet  to  be  successfully  applied  to  boundary  layers 
in  forward  transition.  In  regard  to  the  forced  transition  k-e  calculation  this 
was  run  with  the  k-e  equation  solved  throughout  the  entire  flow  field.  The 
results  indicate  that  for  compressor  cascade  flow  fields  a  model  capable  of 
predicting  transition  accurately  may  be  required  to  obtain  accurate  suction 
surface  velocity  profile  predictions.  Therefore,  based  upon  these  results  two 
possible  approaches  for  further  turbulence  model  work  appear  viable.  The  first 
would  be  to  concentrate  upon  application  of  the  k-e  model  to  flows  in  forward 
transition.  The  second  would  investigate  alternate  k-e  type  models  such  as 
those  discussed  in  Ref.  29. 


Unsteady  Boundary  Conditions 

To  date,  the  Navier-Stokes  cascade  analysis  has  focused  upon  steady  or 
unsteady  flows  through  a  cascade  of  stationary  airfoils.  Another  important 
problem  area  is  that  of  flow  through  a  cascade  in  which  the  blades  are  moving 
relative  one  to  another.  This  problem  arises  in  connection  with  free-vibration 
or  flutter  phenomena  in  the  compressor  or  turbine  stages  of  modern 
turbomachinery  or  engine  configurations.  Early  approaches  to  this  problem  were 
based  upon  classical  linear  analyses  which  assumed  flat  plate  blades  at  zero 
incidence  to  the  mean  flow  (e.g.  Ref.  33).  Although  these  procedures  can  give 
useful  insight  into  the  fluid  dynamic  process,  they  are  obviously  limited  by 
their  assumptions  of  zero  mean  incidence,  zero  thickness  and  camber  and 
inviscid  flow.  More  recent  analyses  although  confined  to  inviscid  flows  have 
removed  the  zero  mean  incidence  and  flat  geometry  restrictions  (e.g.  Ref.  34). 
Although  these  still  treat  the  time-dependent  flow  as  a  linear  perturbation 

32.  Jones,  W.P.  and  Launder,  B.E.:  Some  Properties  of  Sink  Flow  Turbulent 
Boundary  Layers,  Journal  of  Fluid  Mechanics,  Vol.  46,  1972,  pp.  337-351. 

33.  Whitehed,  D.S.:  Vibration  and  Sound  Generation  in  a  Cascade  of  Flat 
Plates  in  Subsonic  Flow.  Aeronautical  Research  Council  R  and  M  3685, 

1970. 

34.  Verdon,  J.M.  and  Caspar,  J.R.:  Development  of  an  Unsteady  Aerodynamic 
Analysis  for  Finite-Reflection  Subsonic  Cascades.  NASA  CR-3455,  1981. 
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from  the  mean  flow,  they  represent  a  more  general  analysis. 

Under  the  present  effort,  initial  consideration  was  given  to  application 
of  the  Navier-Stokes  cascade  analysis  to  the  oscillating  cascade  problem  so  as 
to  reduce  the  inviscid  and  small  perturbation  assumptions.  Since  the  present 
procedure  solves  the  time-dependent  equations  in  a  general  coordinate  system 
which  may  be  a  function  of  time,  the  procedure  is  clearly  applicable  to  the 
oscillating  cascade  problem.  However,  boundary  condition  specification 
presents  a  difficulty. 

One  approach  in  principle  which  is  not  practical  would  consider  all  N 
blades  of  an  N-bladed  cascade  with  the  azimuthal  coordinate  of  the  cascade 
varying  between  0°  and  360°.  Since  the  points  at  0°  and  360°  are  physically 
the  same  point,  a  true  periodic  boundary  condition  could  be  applied. 

Alternately  an  N-bladed  recti-linear  cascade  could  be  considered.  If  N  were 
large,  then  the  flow  in  the  middle  passage  might  be  relatively  insensitive  to 
the  boundary  conditions  applied  at  stagnation  lines  and  wake  lines  for  blades  1 
and  N.  Unfortunately,  either  of  these  options  would  be  impractical  due  to  the 
large  number  of  grid  points  required  and,  therefore,  an  alternative  formulation 
is  sought  whereby  the  flow  in  a  given  passage  or  between  two  adjacent  flow 
periodic  lines  in  the  steady  case  represents  the  flow  in  any  passage  (or 
between  any  two  adjacent  flow  periodic  lines)  in  the  cascade.  Obviously,  if 
the  flow  field  is  to  exhibit  such  periodicity,  the  driving  blade  motion  must 
exhibit  a  periodicity.  The  general  type  of  blade  motion  can  be  given  for  the 
m*1*1  blade  by 


xd  =  ^xfx^i  t  +  <£|  +  rn(cr) 

yd  =  Vy(cV  +  (^2  +  m2<T) 


where  x^  and  y<j  are  displacement  of  the  blade  from  its  mean  position  a  is  a 
constant  interblade  phase  angle,  and  f  represents  a  periodic  function  and  <5X 
and  <5y  represent  oscillation  amplitudes.  This  represents  a  moving  wave 
through  the  blade  row. 
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Turning  to  Fig.  1,  the  problem  can  be  illustrated  by  the  question  of  what 
boundary  conditions  on  boundaries  BA  and  CD  will  represent  the  motion  specified 
by  Eq.  (16)  of  all  blades  external  to  this  portion  of  the  flow  field.  Since  a 
periodic  solution  with  phase  lag  is  sought,  then  the  dependent  variables  on 
corresponding  points  of  lines  AB  and  CD  must  bear  a  periodic  type  relation.  If 
the  value  of  any  dependent  variable,  ¥,  obtained  for  a  steady  solution  is 
denoted  and 


(^)  =  y-xj, 


(17) 


Then  since  these  represent  corresponding  periodic  points.  For  the 

unsteady  case,  the  solution  sought  is  one  in  which  the  disturbed  variable, 
(x|>l)  is  periodic  with  a  possible  interblade  phase  lag  between  corresponding 
points;  i.e.  if 


"Vo  =  ,(",  +  V 


(18) 


then 


WVl  =  *  +  </>[_)  =  f(wt  +  4>g  +  cr) 


The  problem  obviously  arises  as  to  the  choice  of  the  periodic  function  f 
and  the  phase  angle  <|>q  which  may  be  a  function  of  the  spatial  coordinates. 

These  must  be  determined  from  some  other  source  as  they  are  boundary  conditions 
to  the  viscous  problem.  One  possible  method  would  determine  f  and  <(>q  from  a 
small  perterbation  solution  such  as  that  of  Ref.  34.  For  example,  if  the  blade 
motion  for  blade  m  is  given  by 


xd 


=  0 


yd  =  cos(cut  +  me) 


(19) 
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where  £]_,  u)  and  a  are  specified  and  the  inviscid,  small  perturbation  solution 
for  the  dependent  variable  ¥  at  point  G  is 

¥  =  \fr  +  (v//,)G  =  \f/  +  cos(a)t  +  <£Q)  (2°) 

then  at  point  L 

^  ^  +  (v//()L  =  vj/  +  €2  cos(cot  +  -  or )  (21) 

and  a  periodic  type  boundary  condition  of  the  type 

“  v//j  coslcut  +  <£g  -cr)  =  ^  -  i// j  cos((ut  +  <jbQ  )  (22) 

results.  If  <J>q  is  specified  for  each  point  on  the  periodic  boundary,  then 
Eq,  (22)  represents  a  possible  boundary  condition  relating  variables  G  and  L, 
The  existing  Navier-Stokes  cascade  code  has  been  extended  to  allow  a  boundary 
condition  of  the  form  of  Eq.  (22). 

A  second  possible  approach  follows  that  of  Erdos,  Alzner  and  Kalben  in 
their  application  of  the  Euler  equations  to  cascade  problem  (Ref.  35).  In  this 
approach,  the  equations  are  solved  in  the  passage  region  between  points  Lf  and 
Gf  where  L1  and  Gff  are  periodic  points  in  steady  flow  as  are  G*  and  Lf f  (see 
Fig.  1).  In  this  approach,  the  dependent  variables  at  G1  are  related  to  those 
at  Lf  r  at  some  earlier  time  depending  upon  interblade  phase  angle.  Similarly, 
the  dependent  variables  at  Lf  are  related  to  those  at  G,?.  With  this  approach, 
the  second  sweep  equations  which  solve  along  lines  F  G1  and  M  Lf  are  solved 
with  specified  function  boundary  conditions.  This  approach  will  require  some 
run  time  to  establish  the  moving  wave  disturbance  type  flow.  In  addition, 
application  of  these  boundary  conditions  may  adversely  affect  stability. 
However,  the  formulation  does  represent  a  possible  approach. 

CONCLUDING  REMARKS 

The  present  report  represents  continuation  of  an  effort  aimed  at  applying 
the  Navier-Stokes  equations  to  subsonic  and  transonic,  turbulent  cascade  flow 
fields.  During  previous  efforts,  work  focused  upon  development  of  a  cascade 
computational  coordinate  system,  calculation  of  laminar  and  turbulent 


35.  Erdos,  J.,  Alzner,  A.  and  Kalben,  P.:  Computation  of  Steady  and  Periodic 
Two-Dimensional  Non-linear  Transonic  Flow  Problems  in  Turbomachinery, 
(T.C.  Adamson  and  M.F.  Platzer,  Editors),  Hemisphere  Publishing, 
Washington,  1977. 
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cascade  flows,  development  of  a  shock  capturing  technique  and  comparison  of 
calculated  flow  fields  with  data.  Under  the  present  effort,  the  computer  code 
was  revised  to  reduce  run  time,  a  turbulence  energy-turbulence  dissipation 
model  was  added  to  the  code,  a  calculation  was  run  to  convergence  with  this  new 
turbulence  model  and  consideration  given  to  the  forced  blade  vibration 
problem.  The  revised  code  runs  more  than  twice  as  fast  as  the  original  code 
with  little  loss  of  code  flexibility.  It  is  estimated  that  further 
restructuring  would  reduce  run  time  by  an  additional  factor  of  two  without  loss 
of  generality.  The  calculation  run  with  the  turbulence  energy-turbulence 
dissipation  model  reached  convergence  with  no  problem.  The  predicted  surface 
pressure  distribution  agreed  well  with  data,  however,  discrepancies  were 
present  between  the  calculated  and  measured  suction  surface  boundary  layer 
profiles.  One  possible  source  of  discrepancy  is  the  laminar  to  turbulent 
transition  process  and  future  work  should  focus  upon  this  phenomenon.  Finally, 
consideration  has  been  given  to  the  problem  of  boundary  conditions  for  the 
forced  vibration  problem. 
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Fig.  1  -  Cascade  coordinate  grid  -  all  grid  points  not  included. 
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Fig.  2  -  Boundary  layer  profile  (Coles  and  Hirst, 
1968  Stanford  Conference,  Vol.  2, 
Compiled  Data,  ID-1400,  x  =  2.887m) 
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Fig.  3  -  Wall  static  pressure  distribution  -  normal  shock  in  tube. 
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Fig.  4  -  Velocity  profile  -  normal  shock  in  tube. 
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Pig,  5  -  Airfoil  surface  pressure  distribution  for  Hobbs  Cascade. 
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Fig.  7  -  Pressure  surface  boundary  layer  profile  for  Hobbs  cascade. 
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Fig.  8  -  Suction  surface  boundary  layer  profile  for  Hobbs  cascade. 
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Fig.  9  -  Suction  surface  boundary  layer  profile  for  Hobbs  cascade. 
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APPENDIX  A  -  COORDINATE  SYSTEM  GENERATION 


In  brief,  the  coordinate  system  consists  of  a  set  of  two  families  of 
curves;  the  £  *  constant  curves  such  as  lines  FG  or  HI  in  Fig.  1  and  the 
n  =  constant  curves  such  as  ABCD  or  A*EDf  in  Fig.  1.  The  coordinate  system  is 
constructed  by  first  forming  the  inner  loop  A’EDf  which  includes  the  blade. 

The  blade  may  be  either  specified  by  an  analytic  equation  or  by  discrete  data 
points.  If  an  equation  is  used,  then  construction  of  the  inner  loop  is 
straight-forward.  If  the  blade  is  specified  by  discrete  data  points,  then,  in 
general,  the  points  required  on  the  inner  loop  will  not  coincide  with  any  point 
used  for  blade  specification.  In  this  case,  a  curve  fit  is  used  to  obtain  the 
required  inner  loop  points.  The  curve  fit  used  is  based  upon  a  local  parabolic 
fit.  For  any  given  point  required  on  the  inner  loop,  a  parabola  is  fit  through 
three  adjacent  specifying  points,  two  on  the  right  and  one  on  the  left.  A 
second  parabola  is  then  fit  through  the  two  points  on  the  left  and  one  on  the 
right.  The  location  of  the  required  point  is  obtained  via  a  weighted  average 
of  these  curve  fits  with  the  weighting  factor  being  determined  by  the  distance 
from  the  required  point  to  the  center  specifying  point  of  each  parabola.  This 
is  followed  by  constructing  an  outer  loop  ABCD  which  consists  of  periodic  lines 
AB  and  CD  and  a  frontal  curve  BC.  Both  the  inner  and  outer  loops  are  then 
represented  by  parametric  curves 

X  =  x(s),  y  =  y(s) 

where  the  parameter  varies  from  zero  to  unity.  The  present  coordinate 
generation  process  utilizes  a  multi-part  transformation.  First  x  and  y  are 
expressed  as  a  function  of  s?,  the  physical  distance  along  the  curve.  Then  s1 
is  related  to  s  via  a  hyperbolic  tangent  parametrization  centered  about  the 
leading  edge  for  the  inner  loop,  and  a  cubic  polynomial  representation  for  the 
outer  loop.  The  inner  loop  transformation  parametric  representation  is  chosen 
so  as  to  have  rapid  variation  in  the  airfoil  leading  edge  region.  Both  the 
inner  loop  hyperbolic  tangent  transformation  and  the  outer  loop  cubic 
transformation  are  performed  between  s’/2  and  0  and  between  s’/2  and  s?  on  each 
loop.  This  ensures  that  the  distance  in  transformed  space  s  between  any  point 
on  the  upper  periodic  line,  line  AB,  and  point  A  is  the  same  as  the  distance  in 
transformed  space  between  the  corresponding  point  on  the  lower  periodic  line, 
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line  CD,  and  point  D.  Similarly,  corresponding  points  on  the  branch  cut  will 
be  equi-distant  in  s  from  the  branch  cut  end  points,  points  Af  and  D1.  This 
property  is  required  if  corresponding  branch  cut  points  (and  corresponding 
periodic  line  points)  are  to  fall  on  the  same  pseudo-radial  line. 

If  N  pseudo-radial  lines  are  to  be  used,  the  grid  points  on  each  loop  are 
chosen  at  values 


Sj  =  (i-  l)/(  N-l)  i  =  0,  N-l 

and  points  having  the  same  value  of  Sj_  on  the  inner  and  outer  loops  define 
the  end  points  for  a  given  pseudo-radial  line  (e.g.  HI).  Since  the  loops  are 
parameterized  so  that  s  varies  rapidly  in  the  region  of  the  leading  edge,  this 
process  effectively  packs  points  into  the  leading  edge  region.  These  points 
are  then  connected  via  cubic  curve  polynomials  which  allow  specification  of 
the  slope  at  each  end  point  (Ref.  7,  SRA  Rpt.).  For  pseudo-radial  lines  which 
intersect  the  inner  loop  downstream  of  a  user  specified  location  (usually 
x/c  <  .05),  the  slopes  are  specified  so  as  to  be  vertical  or  normal  at  the 
inner  loop  end  points.  For  pseudo-radial  lines  whose  intersection  with  the 
inner  loop  is  upstream  of  this  location  a  smooth  variation  of  the  angle  to  the 
horizontal  with  distance  around  the  leading  edge  is  specified.  With  this 
process,  the  pseudo-radial  lines  are  vertical  as  they  pass  through  the  periodic 
lines  (BA  and  CD)  and,  therefore,  the  metric  data  will  be  continuous  across 
these  lines.  This  process  allows  construction  of  all  pseudo-radial  lines. 

The  second  set  of  coordinate  lines,  the  pseudo-azimuthal  lines  (e.g.  JPK) 
are  obtained  by  normalizing  the  length  along  each  line  with  a  hyperbolic 
tangent  transformation  that  packs  points  near  the  airfoil  surface  (or  branch 
cut).  Then  if  M  pseudo-azimuthal  lines  are  to  be  used,  the  grid  point  for  the 
jth  pseudo-azimuthal  point  occurs  at  rj  =  (j  —1 ) / (M— 1 )  for  each 
pseudo-radial  line.  The  actual  systems  used  in  the  present  effort  consisted  of 
113  points  along  each  loop  and  30  points  along  each  pseudo-radial  line.  The 
plot  shown  in  Fig.  1  is  a  computer  generated  plot  for  the  Jose  Sanz  cascade  in 
which  all  grid  points  are  not  included.  For  clarity,  radial  lines  in  the 
vicinty  of  the  leading  edge  were  omitted  and  loops  in  the  vicinity  of  the  blade 
surface  were  omitted,  leading  to  the  non-smooth  appearance  of  the  grid. 
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The  actual  grid  spacing  normal  to  the  surface  was  4  x  10"^  chords  at  the 
blade  surface  and  the  spacing  along  the  blade  surface  was  4.6  x  10“^  chords 
at  the  blade  leading  edge.  In  addition,  the  grid  extended  two  chords 
downstream  of  the  airfoil  trailing  edge. 

In  generating  the  coordinate  system,  it  is  usually  necessary  to  adjust  the 
parameters  of  the  inner  and  outer  loop  transformation  to  control  the  location 
of  the  pseudo-radial  lines,  to  maintain  a  smooth  variation  of  metric  data  and 
to  minimize  coordinate  nonorthogonality.  In  some  cases  it  may  be  necessary  to 
add  additional  parameter izations  to  provide  a  suitable  grid.  In  general,  the 
additional  transformations  used  are  local  transformations  which  only  affect 
specific  regions  of  the  loops.  For  example,  if  it  is  desired  to  increase  the 
parametric  variation  between  points  A  and  B,  then 

s‘  =  s  s  <  sA 

s'  =  S  +  -|-[l -CO$(ir{S-SA)/(SB  -SA))]  SA  <  S  <  SB 

s'  =  s  +  €  s  >  sa 


would  be  a  suitable  local  transformation. 
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APPENDIX  B  -  SOLUTION  PROCEDURE  [23] 

Background 

The  solution  procedure  employs  a  consistently-split  linearized  block 
implicit  (LBI)  algorithm  which  has  been  discussed  in  detail  in  [12,  21]. 

There  are  two  important  elements  of  this  method: 

(1)  the  use  of  a  noniterative  formal  time  linearization  to 
produce  a  fully-coupled  linear  multidimensional  scheme  which 
is  written  in  "block  implicit"  form;  and 

(2)  solution  of  this  linearized  coupled  scheme  using  a  consistent 
"splitting"  (ADI  scheme)  patterned  after  the  Douglas-Gunn  [22] 
treatment  of  scalar  ADI  schemes. 

The  method  is  thus  referred  to  as  a  split  linearized  block  implicit  (LBI) 
scheme.  The  method  has  several  attributes: 

(1)  the  noniterative  liearization  is  efficient; 

(2)  the  fully-coupled  linearized  algorithm  eliminates  instabilities 
and/or  extremely  slow  convergence  rates  often  attributed  to  methods 
which  employ  ad_  hoc  decoupling  and  linearization  assumptions  to 
identify  nonlinear  coefficients  which  are  then  treated  by  lag  and 
update  techniques; 

(3)  the  splitting  or  ADI  technique  produces  an  efficient  algorithm 
which  is  stable  for  large  time  steps  and  also  provides  a  means  for 
convergence  acceleration  for  further  efficiency  in  computing  steady 
solutions ; 

(4)  intermediate  steps  of  the  splitting  are  consistent  with  the 
governing  equations,  and  this  means  that  the  "physical"  boundary 
conditions  can  be  used  for  the  intermediate  solutions.  Other 
splittings  which  are  inconsistent  can  have  several  difficulties  in 
satisfying  physical  boundary  conditions  [21]. 

(5)  the  convergence  rate  and  overall  efficiency  of  the  algorithm  are 
much  less  sensitive  to  mesh  refinement  and  redistribution  than 
algorithms  based  on  explicit  schemes  or  which  employ  ad  hoc 
decoupling  and  linearization  assumptions.  This  is  important  for 
accuracy  and  for  computing  turbulent  flows  with  viscous  sublayer 
resolution;  and 
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(6)  the  method  is  general  and  is  specifically  designed  for  the 

complex  systems  of  equations  which  govern  multiscale  viscous  flow 
in  complicated  geometries. 

This  same  algorithm  was  later  considered  by  Beam  and  Warming  [36],  but  the 
ADI  splitting  was  derived  by  approximate  factorization  instead  of  the 
Douglas-Gunn  procedure.  They  refer  to  the  algorithm  as  a  "delta  form" 
approximate  factorization  scheme.  This  scheme  replaced  an  earlier  non-delta 
form  scheme  [37],  which  has  inconsistent  intermediate  steps. 

Spatial  Differencing  and  Artificial  Dissipation 

The  spatial  differencing  procedures  used  are  a  straightforward  adaption 
of  those  used  in  [12]  and  elsewhere.  Three-point  central  differene  formulas 
are  used  for  spatial  derivatives,  including  the  first-derivative  convection 
and  pressure  gradient  terms.  This  has  an  advantage  over  one-sided  formulas 
in  flow  calculations  subject  to  "two  point"  boundary  conditions  (virtually 
all  viscous  or  subsonic  flows),  in  that  all  boundary  conditions  enter  the 
algorithm  implicitly.  In  practical  flow  calculations,  artificial  dissipation 
is  usually  needed  and  is  added  to  control  high-frequency  numerical 
oscillations  which  otherwise  occur  with  the  central-dif f erence  formula. 

In  the  present  investigation,  artificial  (anisotropic)  dissipation  terms 
of  the  form 

,  a2 

d  3  u 

2  J[_  _ k  (1) 

j  h.2  3x. 2 

3  3 

are  added  to  the  right-hand  side  of  each  (k-th)  component  of  the  momentum 
equation,  where  for  each  coordinate  direction  x j  ,  the  artificial 
diffusivity  dj  is  positive  and  is  chosen  as  the  larger  of  zero  and  the 
local  quantity  ye  (a  Re^x“l)/Re.  Here,  the  local  cell  Reynolds  number 
Re^xj  for  the  j-th  direction  is  defined  by 

ReAxj  =  Re  |  puj  j  Axj  /  ue  (2) 

This  treatment  lowers  the  formal  accuracy  to  0  (Ax),  but  the  functional  form 
is  such  that  accuracy  in  representing  physical  shear  stresses  in  thin  shear 
layers  with  small  normal  velocity  is  not  seriously  degraded.  This  latter 
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property  follows  from  the  anisotropic  form  of  the  dissipation  and  the 
combination  of  both  small  normal  velocity  and  small  grid  spacing  in  thin 
shear  layers . 

Split  LBI  Algorithm 
Linearization  and  Time  Differencing 

The  system  of  governing  equations  to  be  solved  consists  of  three/four 
equations:  continuity  and  two/three  components  of  momentum  equation  in 

three/four  dependent  variables:  p,  u,  v,  w.  Using  notation  similar  to  that 
in  [12],  at  a  single  grid  point  this  system  of  equations  can  be  written  in 
the  following  form: 

3H(*)/3t  =  D(  <j>)  +  S(  (j>)  (3) 

where  <{>  is  the  column-vector  of  dependent  variables,  H  and  S  are  column- 
vector  algebraic  functions  of  <f>,  and  D  is  a  column  vector  whose  elements  are 
the  spatial  differential  operators  which  generate  all  spatial  derivatives 
appearing  in  the  governing  equation  associated  with  that  element. 

The  solution  procedure  is  based  on  the  following  two-level  implicit 
time-difference  approximations  of  (3): 

(Hn+1-  Hn)/At  =  3(Dn+1+  Sn+1)  (1-6)  (Dn  +  Sn)  (4) 

where,  for  example,  Hn+1  denotes  H(<{>n+1)  and  At  =  tn+*  -  tn.  The 

parameter  3  (0.5-3-  1)  permits  a  variable  time-centering  of  the  scheme, 
with  a  truncation  error  of  order  [At^,  (3  -  1/2)  At]. 

A  local  time  linearization  (Taylor  expansion  about  <)>n)  of  requisite 
formal  accuracy  is  introduced,  and  this  serves  to  define  a  linear  differen¬ 
tial  operator  L  (cf.  [12])  such  that 

D  =  D  +  L  ($  -  $  )  +  0(  At  )  (5) 

Similarly, 

Hn+i  =  Hn+  ( 3H/3cf>)n  (<}>n+1  -  <f>n)  +  0  (At2)  (6) 

sn+i  =  gn+  (9s/3()))n  ((})n+1  -  4>n)  +  0  (At2)  (7) 

Eqs .  (5-7)  are  inserted  into  Eq.  (4)  to  obtain  the  following  system  which  is 
linear  in  <j>n+l 

(A  -  gAt  Ln)  (<)>n+1  -  <f>n)  =  At  (Dn  +  Sn)  (8) 
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and  which  is  termed  a  linearized  block  implilcit  (LBI)  scheme.  Here,  A 
denotes  a  matrix  defined  by 

A  =  ( 3H/34>)n  -  0At  ( 3S/ 3<f>)n  (9) 

Eq.  (8)  has  0  (At)  accuracy  unless  H  =  <|>,  in  which  case  the  accuracy  is  the 
same  as  Eq.  (4). 

Special  Treatment  of  Diffusive  Terms 

The  time  differencing  of  diffusive  terms  is  modified  to  accomodate 
cross-derivative  terms  and  also  turbulent  viscosity  and  artificial  dissipa¬ 
tion  coefficients  which  depend  on  the  solution  variables.  Although  formal 
linearization  of  the  convection  and  pressure  gradient  terms  and  the  resulting 
implicit  coupling  of  variables  is  critical  to  the  stability  and  rapid  con¬ 
vergence  of  the  algorithm,  this  does  not  appear  to  be  important  for  the 
turbulent  viscosity  and  artificial  dissipation  coefficients.  Since  the 
relationship  between  ye  and  dj  and  the  mean  flow  variables  is  not  conven¬ 
iently  linearized,  these  diffusive  coefficients  are  evaluated  explicitly  at 
tn  during  each  time  step.  Notationally ,  this  is  equivalent  to  neglecting 
terms  proportional  to  3  ye/ 3<f>  or  3dj/3<j>  in  Ln,  which  are  formally  pre¬ 
sent  in  the  Taylor  expansion  (5),  but  retaining  all  terms  proportional  to 
ye  or  dj  in  both  Ln  and  Dn. 

It  has  been  found  through  extensive  experience  that  this  has  little  if 
any  effect  on  the  performance  of  the  algorithm.  This  treatment  also  has  the 
added  benefit  that  the  turbulence  model  equations  can  be  decoupled  from  the 
system  of  mean  flow  equations  by  an  appropriate  matrix  partitioning  (cf. 

[21])  and  solved  separately  in  each  step  of  the  ADI  solution  procedure.  This 
reduces  the  block  size  of  the  block  tridiagonal  systems  which  must  be  solved 
in  each  step  and  thus  reduces  the  computational  labor • 

In  addition,  the  viscous  terms  in  the  present  formulation  include  a 
number  of  spatial  cross-derivative  terms.  Although  it  is  possible  to  treat 
cross-derivative  terms  implicitly  within  the  ADI  treatment  which  follows,  it 
is  not  at  all  convenient  to  do  so;  and  consequently,  all  cross-derivative 
terms  are  evaluated  explicitly  at  tn.  For  a  scalar  model  equation 
representing  combined  convection  and  diffusion,  it  has  been  shown  by  Beam  and 
Warming  [38]  that  the  explicit  treatment  of  cross-derivative  terms  does  not 
degrade  the  unconditional  stability  of  the  present  algorithm.  To  preserve 
notational  simplicity,  it  is  understood  that  all  cross— derivative  terms 
appearing  in  Ln  are  neglected  but  are  retained  in  Dn.  It  is  important  to 
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note  that  neglecting  terms  in  Ln  has  no  effect  on  steady  solutions  of  Eq.  (8), 
since  <J>n+l  -  <|>n  =  0,  and  thus  Eq.  (8)  reduces  to  the  steady  form  of  the 
equations:  Dn  +  Sn  =  0.  Aside  from  stability  considerations,  the  only 

effort  of  neglecting  terms  in  Ln  is  to  introduce  an  0  (At)  truncation  error. 

Consistent  Splitting  of  the  LBI  Scheme 

To  obtain  an  efficient  algorithm,  the  linearized  system  (8)  is  split  using 
ADI  techniques.  To  obtain  the  split  scheme,  the  multidimensional  operator  L  is 
rewritten  as  the  sum  of  three  "one-dimensional”  sub-operators  (i  =  1,  2,  3) 
each  of  which  contains  all  terms  having  derivatives  with  respect  to  the  i-th 
coordinate.  The  split  form  of  Eq.  (8)  can  be  derived  either  as  in  [12,  21]  by 
following  the  procedure  described  by  Douglas  and  Gunn  [22]  in  their 
generalization  and  unification  of  scalar  ADI  schemes,  or  using  approximate 
factorization.  For  the  present  system  of  equations,  the  split  algorithm  is  given 
by 

(A  -  gAtL®)  (4>*  -  <]>n)  =  At  (Dn  +  Sn)  '  (10a) 

(A  -  gAtL*)  (♦**  -  4>n)  =  A  (<f>*  -  <j>n)  (10b) 

(A  -  gAtL^)  Un+1  -  4>n)  =  A  (<]>**  -  <J>n)  (10c) 

where  <J>*  and  $**  are  consistent  intermediate  solutions.  If  spatial 
derivatives  appearing  in  and  D  are  replace  by  three-point  difference 
formulas,  as  indicated  previously,  then  each  step  in  Eqs.  (lOa-c)  can  be  solved 
by  a  block-tridiagonal  elimination. 

Combining  Eqs .  ( lOa-c )  gives 

(A  -  BAtL®)  A-1  (A  -  gAtL^)  A_1  (A  -  gAtL^)  ( <f>n+1  -  <t>n)  -  At  (Dn  +  Sn)  (11) 

which  approximates  the  unsplit  scheme  (8)  to  0  (At^).  Since  the  intermediate 
steps  are  also  consistent  approximations  for  Eq.  (8),  physical  boundary 
conditions  can  be  used  for  <J>*  and  <)>**  [12,  21].  Finally,  since  the  L-j_ 
are  homogeneous  operators,  it  follows  from  Eqs.  (lOa-c)  that  steady  solutions 
have  the  property  that  $n+l  =<)>*  =  $**  =  and  satisfy 

Dn  +  Sn  -  0  (12) 

The  steady  solution  thus  depends  only  on  the  spatial  difference  approximations 
used  for  (12),  and  does  not  depend  on  the  solution  algorithm  itself. 
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